function [vol_gy,m1cond_y,m2cond_y,EFGX] = moment_quarterly_output_growth_vol_Nxz(parms,gesol,Phi_Nxz,prob_Nxz)
% see notes
% prob_Nz is the stationary probability

    logfun = @(x)log(x);
    logfun2 = @(x)(log(x)).^2;
    
    Y_NLz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]).*permute(repmat(exp(parms.grid_z(:)),[1 parms.NN parms.Nx]),[2 3 1]);
    
    % 1. compute E[log(X(t)+X(t+1)+X(t+2))|N(t),H(t),z(t)]
    m1cond_y = moment_EFX_conditional_Model20220926(parms,gesol,logfun,Y_NLz);
    
    % 2. compute E[log(X(t)+X(t+1)+X(t+2))^2|N(t),H(t),z(t)]
    m2cond_y = moment_EFX_conditional_Model20220926(parms,gesol,logfun2,Y_NLz);
    
    % 3. compute E[log(X(t)+X(t+1)+X(t+2))*log(X(t+3)+X(t+4)+X(t+5))|N(t),H(t),z(t)]
    EFGX = moment_EFGX_conditional_Model20220926(parms,gesol,@(x)log(x),Y_NLz,reshape(Phi_Nxz*m1cond_y(:),[parms.NN parms.Nx parms.Nz]));
    
    vol_gy = sqrt(max(2*sum(prob_Nxz(:).*m2cond_y(:)) - 2*sum(prob_Nxz(:).*EFGX(:)),0));

end